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Measurements  from  angular  rate  sensors  and  accelerometers  allow 
inertial  navigation.  While  this  report  focuses  on  accelerometers 
(since  they  are  conceptually  simpler),  all  the  discussion  is  directly 
applicable  to  angular  rate  sensors.  In  fact,  estimation  techniques 
were  developed  to  analyze  angular  rate  sensors  attached  to  a  body 
subjected  to  sinusoidal  motion.  System  developers  need  to  be  able 
to  predict  performance  of  a  system,  based  on  the  component 
specifications  supplied  by  the  manufacturers  and  signal  processing. 
Given  knowledge  of  the  noise  level,  the  time  correlation  of  the 
noise,  the  number  of  wavelengths  processed,  and  the  sampling  rate, 
a  model  is  proposed  to  predict  the  accuracy  of  the  estimated 
parameters  of  a  sinusoidal  input.  It  is  hoped  that  this  model  will 
allow  the  designer  to  quickly  determine  the  required  specifications 
for  a  particular  application. 
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FREQUENCY  ESTIMATION  FROM  ACCELEROMETER  MEASUREMENTS 


1.  Introduction 


Measurements  from  angular  rate  sensors  and  accelerometers  allow  inertial 
navigation.  While  this  report  focuses  on  accelerometers  (since  they  are 
conceptually  simpler),  all  the  discussion  is  directly  applicable  to  angular  rate 
sensors.  In  fact,  estimation  techniques  were  developed  to  analyze  angular  rate 
sensors  attached  to  a  body  subjected  to  sinusoidal  motion.  System  developers 
need  to  be  able  to  predict  the  performance  of  a  system,  based  on  the  component 
specifications  supplied  by  the  manufacturers  and  signal  processing.  Given 
knowledge  of  the  noise  level,  the  time  correlation  of  the  noise,  the  number  of 
wavelengths  processed,  and  the  sampling  rate,  a  model  is  proposed  to  predict  the 
accuracy  of  the  estimated  parameters  of  a  sinusoidal  input.  It  is  hoped  that  this 
model  will  allow  the  designer  to  quickly  determine  the  required  specifications  for 
a  particular  application. 

A  simulation  was  developed  to  create  data  for  the  desired  situations.  Two 
estimates  are  pursued.  First,  an  estimation  technique  that  identifies  the  dominant 
frequency1  in  the  signal  is  discussed.  The  second  frequency  estimation  technique 
introduces  new  estimation  criteria  based  on  an  inner  product.  A  model  of 
estimation  error  is  presented,  and  some  statistical  concerns  are  addressed. 

The  frequency  of  the  sinusoid,  the  power/amplitude,  and  the  phase  are  the 
parameters  of  interest.  Notice  that  when  the  signal  of  interest  is  known  to  be 
sinusoidal,  the  bias  or  offset  from  zero  is  of  no  concern.  Once  the  frequency  is 
found,  the  bias  is  the  average  of  the  signal  over  a  whole  number  of  wavelengths. 
It  does  not  affect  the  estimates  of  the  sinusoidal  parameters.  Any  time  an 
accelerometer  has  a  sinusoidal  input,  the  bias  can  be  estimated.  It  may  be 
possible  to  estimate  bias  for  accelerometers  installed  in  artillery  rounds  under  the 
assumption  that  the  coning  motion  is  sinusoidal.  First,  an  accelerometer  must  be 
aligned  so  that  it  is  in  the  coning  plane,  and  the  bias  must  be  calculated  based  on 
the  assumption  that  the  signal  is  sinusoidal.  Then,  in-flight  calibration  can  be 
achieved  via  a  mechanism  to  rotate  the  accelerometer  into  the  axis  of  motion. 

The  frequency  is  the  parameter  of  greatest  interest.  If  the  estimate  of  frequency  is 
wrong,  our  predictions  will  deviate  at  the  rate  of  the  frequency  estimation  error. 
By  having  a  “bound”  on  the  error  in  frequency  estimation,  it  is  possible  to 
determine  the  duration  of  useful  prediction  for  an  estimate.  Errors  in  phase  will 
cause  a  shift  in  the  sinusoid  and  do  not  change  with  time.  Similarly,  estimates  of 


'The  dominant  frequency  is  the  frequency  that  yields  the  largest  amplitude  or  power  when  the 
signal  is  projected  onto  it. 
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the  power  or  magnitude  of  the  sine  wave  will  be  similar  to  the  actual  wave  and 
will  not  change  over  time. 


2.  Simulation 


The  signal  simulation  is  in  SIMULINK®,  an  icon-based  software  product  that 
allows  system  analysis.  An  accelerometer  model  was  developed  which  consisted 
of  the  input,  a  bias  added  to  the  input,  and  a  time-correlated  noise  process  added 
to  the  input.  A  SIMULINK®  block  for  a  time-correlated  noise  process  was 
developed  for  the  accelerometer  model.  This  simple  model  is  justified  if  the 
accelerometer  is  assumed  to  have  been  calibrated.  Scale  factor  corrections  and 
predetermined  biases  are  assumed  to  be  properly  accounted.  The  model 
represents  all  the  un-modeled  biases  as  a  single  bias  and  adds  noise  to  the  signal. 
A  correlated  noise  process  is  used.  Noise  is  correlated  if  the  current  noise  value 
depends  on  the  value  in  the  recent  past.  The  time  correlation  indicates  the 
duration  of  the  dependence.  A  first  order  differential  equation  can  be  used  to 
model  correlated  noise. 

The  input  to  the  accelerometer  was  chosen  to  be  a  sinusoidal  function.  By  setting 
the  parameters  of  the  simulation  noise  level,  time  correlation,  and  step  size,  one 
can  generate  the  desired  data  sets.  Step  size  will  determine  the  number  of  points 
generated  per  unit  time.  This  simulation  did  not  represent  the  internal  workings 
of  the  accelerometer.  The  model  is  shown  in  Figure  1,  and  Figure  2  shows  the 
accelerometer  block.  Output  from  the  model  includes  the  acceleration,  velocity, 
and  distance.  The  model  blocks  labeled  1/s  are  integrators. 


Distance 


Figure  1.  Signal  Generation  Model. 
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Figure  2.  Accelerometer  Model. 


3.  Signal  Processing 


In  this  case,  both  the  power  spectrum  and  a  Fourier  analysis  are  not  that  helpful 
since  the  quantities  of  interest  are  the  parameters  of  the  sinusoid  with  the  greatest 
power.  In  Fourier  analysis,  an  orthogonal  basis  of  sinusoids  is  used  to  represent 
the  signal.  This  set  will  capture  all  the  energy  of  the  signal;  however,  it  will 
usually  miss  the  actual  frequency  of  the  signal.  The  set  of  orthogonal  sinusoids  is 
discrete  and  if  the  frequency  of  interest  falls  in  between  the  Fourier  frequencies,  it 
is  difficult  to  estimate  precisely. 

A  set  of  routines  was  developed  to  find  the  frequency  of  maximum  power.  The 
method  used  is  similar  to  either  an  ideal  receiver  or  the  method  used  for  wavelets. 
Wavelets  use  an  orthogonal  set  of  functions  to  represent  a  signal.  They  are  useful 
when  a  few  of  the  wavelets  can  capture  most  of  the  information  in  the  signal. 
Generating  a  good  set  of  orthogonal  functions  for  a  particular  problem  can  be  a 
difficult  task.  In  this  case,  the  only  information  of  interest  is  that  associated  with 
the  sinusoid  of  the  greatest  power;  thus,  all  other  information  in  the  signal  can  be 
ignored.  Since  the  signal  is  known  to  be  sinusoidal,  the  time  delay  and  amplitude 
can  be  ascertained  by  projections  of  the  signal  onto  sine  and  cosine  waves.  Note 
that  different  frequencies  require  different  time  intervals  in  order  to  maintain  the 
orthogonality  of  the  sine  and  cosine  functions.  First,  a  routine  to  find  the  power 
and  phase  of  a  specified  frequency  was  developed.  This  routine  was  then  used  to 
map  the  power  surface  within  a  region  of  the  frequency  space.  It  was  observed 
that  power,  as  a  function  of  frequency,  was  a  smooth  function  in  the  region  of 
maximum  power.  An  iterative  search  routine  using  a  quadratic  estimator  was 
designed  to  find  the  frequency  of  maximum  power. 
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A  second  estimation  method  for  frequency  was  also  devised.  This  method  is 
based  on  the  observation  that  sine  and  cosine  functions  of  the  desired  frequency 
are  orthogonal.  The  inner  product  of  the  signal  projected  onto  the  sine  and  cosine 
will  be  zero  if  the  proper  frequency  is  selected.  By  calculating  these  quantities  for 
different  frequencies,  one  can  ascertain  the  desired  frequency.  The  estimation 
process  then  becomes  a  zero  crossing  problem.  This  estimation  criterion  will  be 
referred  to  as  the  “inner  product  criterion.”  As  in  the  previous  case,  only 
frequencies  completing  a  whole  number  of  cycles  over  a  selected  data  set  could 
be  directly  evaluated. 

Methods  to  select  the  parameter  values  can  be  based  on  other  criteria.  Usually,  a 
criterion  that  is  a  function  of  the  residuals  is  selected.  The  parameters  can  be 
chosen  so  that  the  sum  of  the  residuals  is  zero.  Minimization  of  the  HI,  H2,  or  H- 
infinity  norms  of  the  residuals  can  be  achieved.  These  norms  correspond  to  the 
sum  of  the  absolute  value  of  the  residuals,  the  sum  of  the  squares  of  the  residuals, 
and  the  maximum  residual. 


4.  Performance  Model 


In  the  model,  estimation  error  is  assumed  to  be  a  function  of  noise  level,  time 
correlation,  number  of  wavelengths  processed,  and  the  number  of  data  points  per 
wavelength.  The  error  is  expected  to  be  infinite  if  no  datum  is  processed;  thus,  at 
levels  of  zero  for  both  wavelengths  processed  and  number  of  data  points,  the 
predicted  error  should  be  infinite.  Error  is  expected  to  be  zero  if  there  is  no  noise 
in  the  accelerometer;  thus  at  zero  noise  level,  the  estimation  error  should  be  zero. 
A  multiplicative  model  with  unknown  exponents  was  considered  to  be  of  the 
correct  form.  Negative  exponents  will  indicate  that  the  error  becomes  smaller 
when  the  value  of  an  independent  variable  increases,  while  positive  exponents 
will  cause  the  error  to  increase  as  the  magnitude  of  an  independent  variable 
increases.  The  number  of  wavelengths  and  samples  per  wavelength  should  have 
negative  exponents,  while  the  noise  level  should  have  a  positive  exponent.  These 
assumptions  can  be  used  as  a  logical  verification  of  the  model  after  the  estimation 
process.  For  this  case,  least  squares  can  be  used  to  estimate  the  exponents  of  the 
independent  variables  via  the  log  transformation  of  the  estimation  error  (y)  as  the 
dependent  variable.  Using  Xs  to  represent  the  independent  variables,  the  model 

a  p  8  <p 

y  =  cXx  Xi  XiX4s 

transforms  to  the  model 

ln(  y)  =  ln(c)  +  a  ln(^)  +  p  ln(^2)  +  8  ln(j3)  +  <p  ln(_^4)  +  ln(^)  . 

Note  for  this  model  that  if  the  dependent  variable  were  raised  to  a  power,  this 
would  result  in  a  multiplication  of  each  of  the  estimated  parameters.  Typically, 
for  a  statistical  interpretation  to  be  appropriate,  the  residual  error  term  should 
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have  a  Gaussian  distribution.  Statistically,  the  goal  is  to  estimate  the  average 
error  and  be  able  to  make  a  statement  of  the  form  “95%  of  the  errors  are  less  than 
X  for  these  values  of  the  independent  variables.” 


5.  Statistical  Issues 


The  model  will  be  used  in  the  original  domain.  In  the  log  domain,  small  errors 
have  large  (negative)  magnitudes,  while  large  errors  have  smaller  magnitudes.  A 
least  squares  estimation  will  tend  to  give  more  influence  to  numbers  of  greater 
magnitude;  however,  it  is  desirable  for  the  model  to  more  accurately  predict  large 
errors.  When  the  log  model  is  transformed  into  the  original  domain,  the  sum  of 
the  residuals  is  no  longer  zero.  This  can  be  corrected  by  adjusting  the  parameters 
using  a  Newton-Raphson  method  to  find  the  proper  value  of  the  parameters. 
Another  approach  is  to  minimize  a  function  of  the  residuals.  Typically,  the  sum 
of  the  absolute  values,  the  sum  of  the  squares,  or  the  greatest  deviation  is 
considered  as  a  possible  metric  to  be  minimized.  A  collection  of  routines  was 
developed  to  minimize  the  selected  metric. 

After  the  parameters  are  adjusted  to  minimize  a  criterion,  a  chi-square  test  can  be 
used  to  test  the  fit  of  the  model.  Through  the  formation  of  the  ratio  of  the  squared 
estimation  error  over  the  model-predicted  variance,  a  chi-square  random  variable 
is  formed.  By  summing  these  and  comparing  the  sum  to  the  proper  critical 
values,  one  can  test  the  null  hypothesis  (i.e.,  that  the  proposed  model  is  adequate 
for  the  observed  data).  Note  that  the  model  derivation  is  independent  of  the  chi- 
square  test.  Another  possibility  would  be  to  search  for  the  parameters  that  give 
the  best  fit  to  a  chi-square  distribution  (50%  level).  This  latter  method  would  be  a 
form  of  maximum  likelihood  estimation. 


6.  Discussion  of  Estimation  Criteria 


The  inner  product  criterion  and  the  maximum  power  criterion  for  estimating  the 
frequency  were  compared.  Over  the  region  observed,  the  power  curve  had  a 
quadratic  component  (see  Figure  3);  however,  increasing  the  amount  of  data  or 
number  of  wavelengths  processed  also  increased  the  curvature  of  the  power  curve. 
This  increased  curvature  is  not  modeled  well  by  quadratic  functions;  higher  order 
estimators  become  necessary.  As  shown  by  Figure  4,  the  inner  product  plot 
appears  to  be  sinusoidal  and  has  a  steep  slope  when  it  crosses  zero.  The  linear 
model  for  the  inner  product  criterion  remains  valid  as  the  number  of  wavelengths 
processed  increases.  For  the  inner  product  criterion,  an  increase  in  the  number  of 
wavelengths  processed  increases  the  slope  as  it  crosses  zero.  Noise  tends  to 
flatten  the  curves  for  both  criteria;  however,  the  steep  slope  of  the  inner  product 
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7.  Bias 


Even  without  the  addition  of  a  deterministic  bias,  a  bias  is  present  in  the  data. 
There  are  two  elements  that  cause  this  undefined  bias.  First,  the  phase  of  the 
signal  at  the  beginning  of  processing  will  cause  a  small  but  persistent  offset. 
Suppose  the  accelerometer  signal  is  a  sine  wave.  The  integration  of  the  signal 
over  the  first  180  degrees  yields  the  maximum  positive  value,  and  integration  over 
the  next  180  degrees  returns  the  result  to  zero;  thus,  the  velocity  term  is  always 
positive.  See  Figure  5  for  assistance  in  visualizing  the  deterministic  bias.  The 
integral  of  the  velocity  will  then  digress  from  the  true  value.  Through  this 
argument,  it  is  evident  that  the  phase  of  a  sinusoidal  signal  will  cause  a  bias  in  its 
integral.  Next,  recall  that  the  integration  operator  is  known  to  be  linear.  The 
integral  of  a  signal  plus  noise  can  be  evaluated  from  the  sum  of  the  integral  of  the 
signal  (as  previously  discussed)  and  the  integral  of  the  noise.  The  integral  of 
noise  is  a  random  walk.  A  random  walk  has  an  equal  chance  of  increasing  or 
decreasing  at  any  point  in  time;  the  best  estimate  of  its  future  value  is  its  current 
value.  Any  offsets  attributable  to  the  random  walk  cause  a  bias  in  the  velocity 
term.  This  component  of  the  bias  cannot  be  predicted  a  priori  and  will  change  as 
a  function  of  time.  System  performance  will  be  determined  by  how  well  the 
effects  of  this  random  walk  can  be  minimized.  High  levels  of  noise  and  long  time 
correlation  increase  the  magnitude  of  the  bias  caused  by  a  random  walk.  Figures 
6  and  7  provide  examples. 


Figure  5.  No  Noise  Response. 
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Velocity 

Distance 


Figure  7.  High  Noise  Condition 


8.  Discussion  of  Variables 


The  four  variables  considered  (number  of  points  per  wavelength,  number  of 
wavelengths,  noise  strength,  and  time  constant  of  the  noise)  were  investigated  in 
ranges  considered  to  be  useful  for  inertial  navigation. 

Although  it  was  possible  to  estimate  frequency  with  ten  observations  per  cycle, 
the  estimation  methods  were  not  always  stable  in  this  region.  At  20  observations 
per  wavelength,  the  estimation  methods  were  stable,  and  it  was  apparent  that  an 
increase  in  the  number  of  observations  beyond  20  per  wavelength  did  not  improve 
the  estimation  fidelity.  This  led  to  the  conclusion  that  the  variable  (observations 
per  wavelength)  was  asymptotic  and  had  reached  a  near  maximum  value  around 
20  observations  per  wavelength.  In  the  future,  the  estimation  routines  could  be 
enhanced  for  better  performance  in  regions  of  few  observations  per  wavelength. 

The  time  constant  of  the  noise  process  indicates  the  auto-correlation  time  of  the 
noise.  At  high  noise  levels,  the  effect  of  a  long  time  constant  is  to  decompose  the 
sinusoid  into  segments  and  increase  the  bias  in  the  random  walk,  making  it 
difficult  to  estimate  any  of  the  sinusoidal  parameters  precisely.  Based  on 
observations  of  accelerometer  and  angular  rate  sensors,  the  time  correlation  is 
small.  Using  a  time  correlation  of  0.01  or  larger  had  no  discernible  effects  on 
parameter  estimation.  The  model  for  the  time  constant  is  a  first  order  differential 
equation.  In  some  cases,  the  time  constant  can  be  thought  of  as  a  means  of 
removing  the  effects  of  un-modeled  differential  equations. 

The  variables  accounting  for  the  noise  level  and  the  number  of  wavelengths 
processed  were  retained  in  the  model.  The  final  model  containing  these  three 
parameters  is 

a  p 

y  =  cx,  x2£ 

Notice  that  the  noise  (XO  should  have  a  positive  exponent,  and  as  the  noise  value 
goes  to  zero,  the  error  should  also  approach  zero.  The  number  of  wavelengths 
(Xj)  should  have  a  negative  exponent.  The  error  should  approach  infinity  as  the 
number  of  wavelengths  goes  to  zero.  This  provides  for  interesting  behavior 
around  the  point  (0,0).  Fortunately,  the  region  of  interest  does  not  contain  this 
point.  Obviously,  the  model  is  only  considered  for  positive  values  of  noise  and 
number  of  wavelengths  processed.  Figure  8  is  included  to  help  the  reader 
visualize  the  model. 
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Figure  8.  Error  Function. 


9.  Fitting  Model  Parameters 


Using  the  accelerometer  model,  we  collected  60  observations.  Noise  levels  were 
set  in  relation  to  the  sinusoid  amplitude.  At  noise  levels  of  0.2,  0.4,  0.8,  and  1.2, 
the  estimation  routine  processed  two  and  four  wavelengths  of  data.  For  the  0.8 
and  1.2  noise  levels,  six  wavelengths  of  data  were  also  processed.  Six 
observations  were  collected  for  each  pair  of  conditions. 

In  the  first  part  of  the  fitting  process,  the  model  was  transformed  to  the  log 
domain,  and  regression  provided  the  initial  parameter  estimates.  Two  routines 
were  devised  to  search  the  parameter  space.  The  first  routine  calculates  the  HI, 
H2,  H-infinity  norms,  and  the  sum  of  the  residuals.  The  second  routine  calculates 
the  p-value  for  a  chi-square  test,  based  on  the  assumption  that  the  model  gives  the 
standard  deviation  of  the  estimation  process.  With  this  software,  a  heuristic 
search  method  was  used  to  find  a  set  of  parameters  that  were  acceptable  for  both 
the  norm  and  chi-square  criteria. 

The  exponent  of  the  noise  term  was  determined  to  be  close  to  1  and  was  set  to  l. 
The  model  then  is  linear  along  the  noise  dimension.  The  exponent  of  the  variable 
for  wavelength  was  -5/3.  The  decrease  in  error  along  the  wavelength  dimension 
is  dramatic.  The  third  parameter  is  a  factor  that  adjusts  the  overall  height  of  the 
surface;  its  value  is  0.1083. 
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10.  Model  Testing 


Two  assumptions  were  used  to  test  the  model.  First,  estimation  errors  were 
assumed  to  have  a  normal  distribution.  Next,  the  model  was  assumed  to  be  true. 
The  value  given  by  the  model  was  the  standard  deviation  of  the  error  for  that 
point.  Under  these  assumptions,  the  observed  error  over  the  model  standard 
deviation  forms  a  chi-square  variable  with  one  degree  of  freedom  when  it  is 
squared. 

The  model  acceptance  test  consisted  of  20  observations.  The  noise  values  were 
randomly  chosen  from  the  range  0.1  to  1.2;  the  number  of  wavelengths  was  a 
randomly  selected  integer  from  the  interval  (2,7).  The  chi-square  variables  were 
formed  and  summed.  The  ensuing  chi-square  test  with  20  degrees  of  freedom 
gave  a  p-value  of  .38  for  the  model.  The  p-value  is  the  probability,  when  the 
assumptions  are  true,  that  if  the  test  were  repeated,  the  new  test  value  would 
exceed  the  observed  value.  The  observed  value  is  in  the  region  where  we  fail  to 
reject  the  assumptions.  Failure  to  reject  the  assumptions  does  not  imply  that  the 
model  is  correct — only  that  based  on  the  data,  a  rare  event  did  not  take  place.  For 
engineering  purposes,  the  model  was  accepted  as  adequate  for  predicting  the 
standard  deviation  of  the  frequency  estimator. 


11.  Concluding  Remarks 


Accelerometer  bias  was  shown  to  have  both  a  deterministic  and  random 
component.  This  is  in  addition  to  biases  caused  by  electronics,  un-modeled 
temperature  effects,  and  so  forth,  and  indicates  that  despite  all  efforts  to  control  or 
calibrate  an  accelerometer,  there  will  be  a  random  bias  caused  by  the  random 
walk  effect  of  the  integrated  noise.  This  random  walk  effect  is  statistically 
proportional  to  the  strength  of  the  noise.  In  stressful  environments  with  lots  of 
vibration,  it  may  be  a  major  factor  in  accelerometer  bias. 

Although  this  investigation  focused  on  frequency  estimation,  the  same  procedure 
can  be  used  to  estimate  amplitude  and  phase  errors.  These  models  would  be  of 
the  same  form.  Although  the  exact  same  analysis  is  pertinent  for  angular  rate 
sensors,  the  accuracy  of  the  amplitude  estimate  is  extremely  important  in  the 
interpretation  of  angular  rate  measurements. 

The  model  validation  technique  is  an  application  of  statistical  theory  to  quantify 
the  performance  of  a  model.  The  area  of  model  evaluation  is  important  to  all 
areas  of  research.  The  technique  developed  here  should  be  investigated  in  more 
detail.  The  current  use,  while  logically  correct,  reverses  the  traditional  role  of  the 
null  hypothesis. 
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In  this  case,  the  signal  was  a  sinusoid  corrupted  by  noise.  The  approach  taken 
was  to  project  the  signal  onto  sinusoids  and  search  the  parameter  space  to  find  the 
best  set  of  parameters.  Two  methods  were  used  to  define  “best”;  the  first  was  to 
use  maximum  power,  and  the  second  method  used  an  inner  product  criterion.  In 
situations  when  the  researcher  knows  the  functional  form  of  the  signal,  the  signal 
can  be  projected  onto  that  function  and  the  desired  parameters  can  be  defined  as 
those  that  yield  the  most  power.  In  this  case,  the  researcher  can  pretend  that  he  or 
she  has  a  set  of  wavelets,  but  only  the  first  one  will  have  useful  information, 
allowing  the  others  to  be  ignored. 

With  flight  data,  all  frequencies  change  during  the  flight.  Frequency  can  be 
tracked  via  either  of  the  estimation  techniques  discussed.  If  a  model  of  the 
frequency  change  were  available,  a  Kalman  filter  could  be  made  to  increase 
tracking  fidelity.  In  many  applications  requiring  frequency  estimation,  such  as 
radar,  banks  of  filters  are  used,  and  the  filter  with  the  greatest  output  is  assumed 
to  have  the  correct  frequency.  Both  frequency  estimation  techniques  mentioned 
here  can  be  implemented  to  improve  performance  and  reduce  the  number  of  filters 
required. 
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